Critical Casimir forces for Ising films with variable boundary fields 
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Monte Carlo simulations based on an integration scheme for free energy differences is used to 
compute critical Casimir forces for three-dimensional Ising films with various boundary fields. We 
study the scaling behavior of the critical Casimir force, including the scaling variable related to the 
boundary fields. Finite size corrections to scaling are taken into account. We pay special attention to 
that range of surface field strengths within which the force changes from repulsive to attractive upon 
increasing the temperature. Our data are compared with other results available in the literature. 
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I. INTRODUCTION 

Forces induced by thermal fluctuations can be very sen- 
sitive to tiny changes in temperature. This is exemplified 
by effective forces arising between two surfaces confin- 
ing a fluid close to its critical point, for which a slight 
variation in temperature can lead to pronounced changes 
in their range and magnitude. The universal features 
of these so-called critical Casimir forces are captured by 
scaling functions they have been studied theoret- 

ically and experimentally for systems belonging to the 
bulk universality classes of the XY and the Ising model 
The XY model describes quantum fluids, such 
as liquid ^He close to its normal-superfluid phase tran- 
sition or a '^He-^He mixture close to its tricritical point, 
whereas, e.g., a classical binary liquid mixture near its 
demixing point or a simple fluid close to a liquid-gas crit- 
ical point belong to the Ising universality class. 

In the systems studied so far experimentally, the mea- 
sured critical Casimir forces have been either attractive 
or repulsive throughout the whole temperature range. 
(The addition of salt to a critical oil-water mixture 
presents a notable exception in that, under favorable con- 
ditions, on route to the critical demixing point the sign 
of the critical Casimir force can change twice Q due to 
a coupling between the noncritical charge density and 
the critical order parameter field.) Here, we investigate 
simple systems which provide the possibility of changing 
the sign of the critical Casimir forces upon varying the 
temperature. Analytic studies and computer simulations 
indicate that the sign of the critical Casimir force is de- 
termined by the properties of the confining surfaces, i.e., 
by the boundary conditions (BCs) which they impose 
on the fluctuations of the order parameter characterizing 
the underlying second-order phase transition. Indirect 
measurements of the Casimir scalingfunction, inferred 
from wetting films of superfiuids H, Q and of classical 
binary liquid mixtures [1, Q, are consistent with these 
predictions. For pure ^He one has symmetric Dirichlet- 
Dirichlet (O, O) BCs because the quantum mechanical 
wave function of the superfiuid state vanishes at both 



confining interfaces. This gives rise to attractive critical 
Casimir forces [lol - [l^ . For wetting films of ^He-^He mix- 
tures, (-I-, O) BCs are realized because due to quantum 
mechanical effects a ^He-rich layer forms near the solid- 
liquid interface and favors the superfiuid pha se giving 
rise to the so-called surface transition [T^. Il5|: (-I-) in- 
dicates a symmetry-breaking BC with the surface com- 
pletely ordered. Upon reaching the surface transition 
the superfiuid order parameter becomes nonzero at the 
solid surface whereas it vanishes at the fluid-vapor in- 
terface of the wetting film. These asymmetric BCs give 
rise to a repulsive Casimir force. Measurements for wet- 
ting films of certain classical binary liquids mixtures have 
been found to be in agreement with (H — ) BCs corre- 
sponding to a strong opposing preferential adsorption 
of the two s pec ies of the mixture at the two confining 
surfaces [lol. Illl [T6l . Il7| . Within the framework of an 
Ising magnet (which is equivalent to the lattice model of 
a binary mixture) or within the continuum field theory 
for the order parameter, this amounts to the presence 
of strong antagonistic symmetry-breaking surface fields 
Hi and H2 which couple linearly to the order parameter 
and give rise to repulsive critical Casimir forces. Direct 
evidences for critical Casimir forces have been provided 
by studying the Brownian motion of a single colloidal 
particle near a flat substrate surface and immersed in 
the binary liquid mixture of water and lutidine (isl [l9| . 
The experimental results for the cases in which the col- 
loid and the substrate surface preferentially adsorb the 
same species of the mixture are consistent with (H — h) or 

( ) BCs, whereas for cases in which the particle and 

the surface preferentially adsorb different species of the 
mixture the results agree with the occurrence of ( — h) 
or (H — ) BCs. Whereas the theoretical and experimental 
understanding of critical Casimir forces in the presence 
of strong or vanishing surface fields has reached a mature 
level, here we set out to study the infiuence of variable 
weak surface fields. 

Dirichlet and (±) BCs are the renormalization-group 
fixed-point boundary conditions corresponding to the so- 
called ordinary surface universality class (O) and the 
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normal transition surface universality class, respectively 
[20l - |22j . The ordinary transition corresponds to the bulk 
phase transition occurring in the absence of surface fields 
and with a reduced tendency to order at the surface. 
Within a mean field picture the latter is described by a 
surface scaling field c so that 1/c plays the role of an ex- 
trapolation length of the order parameter profile; c = oo 
defines the ordinary transition fixed point. The normal 
transition occurs for systems with strong surface fields 
and which exhibit a reduced tendency to order if these 
surface fields are switched off. The normal transition is 
defined by the fixed point (Hi = oo,c = oo). As indi- 
cated by the nomenclature the normal transition is the 
generic situation for a fluid. In a spin model as discussed 
below Hi = Hi/ J is dimensionless with J as an interac- 
tion constant (see below). 

Near the ordinary transition there is a single linear 
scaling field gi = Hi/ associated with the surface field 
of strength Hi and the surface enhancement parameter 
c The scaling exponent is y = (Af - Af*) /$, 

are the surface counterparts of the bulk 
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where A 

gap exponent A and $ is a crossover exponent [2C1| - |22| . 
For the three-dimensional {D = 3) Ising model one has 
Af ~ 0.46(2) [ilLAf ~ 1.05 [ill, $ ~ 0.68 [HI, and 
?/~0.87;i/~0.63 [2J] is the critical exponent of the bulk 
correlation length ^f, = ,^^|t|~^ with the reduced temper- 
ature t = {13c- l3)/(3 = (T - rc)/Tc; ± corresponds to 
t ^ 0. The corresponding surface scaling variable can 
be chosen as (Co')~^S'i . For i — > the scaling 
variables tend to their fixed point values, and the scal- 
ing functions assume asymptotic forms corresponding to 
the respective fixed points [2l|, (2^. The scaling vari- 
able |(Co')~^5i|i|~^°'^'*r^^°'^'* is proportional to the ratio 
/ £i between the true bulk correlation length and the 
length £i introduced by the scaling field gi : 



(1) 



The fixed point dominated critical regions correspond to 
either the divergence ( (-I-) or (— ) fixed-point BCs) or 
the vanishing ( (O) fixed-point BCs) of this ratio. The 
length £i corresponds to the range of distances from the 
surface within which the order parameter profile responds 
linearly to the presence of a surface field Hi [2^ [2^ . (A 
precise definition of £i will be provided below.) 

Depending on the interplay between £i and the length 
ic = (Cdc)^"/* associated with the surface enhancement 
parameter c one finds various asymptotic regimes for the 
short-distance behavior z <C of the order parameter 
profile [2^ [25I, [13]. At the bulk critical point one has 
4'cri{z) ^ giz'^ for distances £c ^ z <C ^1 from the sur- 
face with K = {A°/'^ - f3)/v and (j)cri[z) - z-f^/" for 
distances £c <C £1 <C z from the surface. We note that 
near the ordinary transition fixed-point (i.e., large c) the 
length Ic is small whereas the length ^l can be large or 
small. Within mean field theory one has k — Q due to 
^ord(-2) = 4) = 1/2 and v{D = 4) = 1/2 whereas one 
has k{D = 3) ~ 0.23 [H and {P/v){D = 3) ~ 0.52 [H]. 



Consequently, the critical order parameter (OP) profile 
turns out to be a nonmonotonic function of z. For z <^ (.1 
the OP increases upon increasing z, at z ~ ii it reaches 
a maximum and only for z 3> ^1 the universal "normal" 
fixed-point behavior, i.e., the decay of the OP occurs [25j . 
Accordingly the position of this maximum can serve as 
a definition for the length €1 [2^, [2^ . With increasing 
surface field strength the surface-near regime with the 
aforementioned increase ^ z'^ of the OP becomes nar- 
rower, and eventually for Hi ^ 00 the length scale £1 
goes to zero, such that this regime disappears and the 
normal transition behavior ^ z^^'^ is attained through- 
out. 

For the 2D Ising model on the square lattice with lat- 
tice constant a, the length £1 has been extracted from 
an exact result for the scaling function of the OP pro- 
file below and above T^, the profile at Tc has not been 
reported. In the case that the exchange coupling be- 
tween spins in the surface row is the same as in the 
bulk, the OP scaling function depends on the scaling 
variable ^b/h with li = (a/2) tanh(X)/(tanh /ii)^, where 
K ~ J/{kBT) is the dimensionless reduced exchange 
couphng between Ising spins and hi = Hi/^ksT) [2^. 
Thus for weak surface fields and in the limit K — > Kc, 
i^{Kc) = {a/2)K-'^tanh{Kc)/Hf = 1.066{A)aHi-'^ = 
1.879(0)^^fff 2^ where = 0.5 ln(l + V2) 2± 0.44 is the 
critical coupling and S,^ = a/(4i^c)- This is in line with 
Eq. © due to iy{D = 2) = 1 and Af'^(L> = 2) = 1/2. 
Examination of the OP profiles for T ^ Tc shows that 
the maximum occurs at Zmax — 1.5^i which implies 
£1 ~ 2.8^0+^1"'. 

Studies of systems belonging to the Ising universality 
class [2§-[3^ showed that near bulk criticality the pres- 
ence of the length scale £1 has important consequences 
for finite-sized systems such as slabs of thickness L. For 
these systems the relevant lengths are the bulk correla- 
tion length the distance L between the two confin- 
ing surfaces which exert fields Hi and H2, and the cor- 
responding lengths £1 and £2- The asymptotic critical 
region, associated with (-I-) or (— ) fixed-point boundary 
conditions at the surfaces i = 1,2, corresponds to L ^ ^i, 
whereas corrections proportional to £i/ L are expected to 
be relevant for L £i. In the crossover regime the critical 
properties of the confined systems are particularly sensi- 
tive to the values of the surface fields, i.e., whether one or 
both length scales £i become comparable to or even larger 
than the distance L, together with L, <C £,b- For exam- 
ple, in films with identical surface fields, i.e.. Hi = H2 
and £1 = £2, at bulk criticality and for weak surface fields 
the order parameter profile exhibits two symmetric max- 
ima at z ~ £1 and z ~ L — £1; for even weaker fields 
so that ^1 ~ L these maxima merge into a single one at 
midpoint z ~ L/2 [2^, [s^. Concomitantly the critical 
Casimir amplitude as a function of the surface field Hi , 
i.e., the critical Casimir force at the bulk critical temper- 
ature, exhibits a maximum absolute value at L ~ £1 [29j 

For symmetric surfaces, the effect of variation of the 
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amplitude of Hi on the temperature dependence of 
the critical Casimir force, i.e., the crossover behavior 
between the ordinary and normal surface universality 
classes, was studied within the two-dimensional (2D) 
Ising model by using the quasi-exact numerical density- 
matrix renormalization-group method [sij and within 
continuum mean- field theory [3^ . For L/£i ~ 1 these 
results show strong deviations of the force scaling func- 
tion from its universal fixed-point behavior such as the 
occurrence of two minima, one above and one below Tc, 
but no change in sign as the temperature is varied. It 
turns out that only strongly asymmetric surface fields 
can lead to, even multiple, sign changes of the criti- 
cal Casimir forces upon varying the temperature. This 
has been demonstrated rigorously for 2D Ising films [1^ , 
within mean- field theory for the same geometry [3^ , and 
it was supported by our preliminary results from Monte 
Carlo simulations of simple cubic Ising slabs [s^] ■ Further 
evidence has been provided by Monte Carlo simulations 
of the improved Blume-Capel model in the film geom- 
etry [s^]. The Blume-Capel model has a second-order 
phase transition which also belongs to the 3-D Ising uni- 
versality class. It offers the opportunity that a careful 
choice of the interaction parameters of this model allows 
one to eliminate leading corrections to finite-size scaling 
(see also Ref. UHl). As it will be discussed below, con- 
trolling finite-size corrections is essential for inferring the 
scaling functions of critical Casimir forces from Monte 
Carlo simulation data. In the following, we shall present 
a Monte Carlo simulation study of the critical Casimir 
forces for the 3-D Ising model in a slab geometry with 
freely variable surface fields applied at its bottom and 
top surfaces. Our scan of the parameter space extends 
the one presented in Ref. HJ. As mentioned above, in 
Ref. [l^l certain preliminary results of this study were 
reported together with a detailed continuum mean-field 
analysis. 

The analytic results and the simulation data for the 
scaling functions of the critical Casimir forces for weak 
surface fields can be probed experimentally and they offer 
application perspectives for soft matter systems such as 
tuning the properties of colloidal suspensions. A first at- 
tempt to investigate experimentally the effects of gradual 
changes in the properties of confining surfaces on criti- 
cal Casimir forces was made recently by studying colloids 
suspended in a critical mixture of water and lutidine j36j . 
These experiments have demonstrated the ability to con- 
tinuously tune the order parameter boundary conditions 
at the confining surfaces. This was achieved by a chemi- 
cal treatment of a solid substrate such that it produces a 
spatial gradient of the adsorption preference for lutidine 
and water molecules. Depending on the position of a sin- 
gle disolved colloidal particle at this structured surface 
a smooth transition from attractive to repulsive critical 
Casimir forces was found. However, these experimental 
observations have not yet been cast into a universal scal- 
ing function of the critical Casimir potential which has 
to change sign as function of the effective surface field. 



Our presentation is organized as follows. In Sec. Ullwe 
introduce our model, define the range of parameters for 
which we perform our computations, and briefly present 
the relevant theoretical background. In Sec. IIIII we de- 
scribe the numerical method employed in order to infer 
the scaling functions of the critical Casimir forces from 
the MC simulation data. In Sec. IIVI we discuss correc- 
tions to scaling which we take into account in order to 
obtain data collapse signalling scaling. Section |V] con- 
tains our results. We provide a summary and conclusions 
in Sec. Ell 



II. MODEL AND THEORETICAL 
BACKGROUND 

In the spirit of the universality of critical phenomena 
we study the simplest representative of the 3-D Ising uni- 
versality class, i.e., the three-dimensional Ising model 
defined on a simple cubic lattice. We consider a slab 
geometry. The dimensionless volume of the system is 
L^x LyX Lz where = Ly ^ and A = x Ly with 
periodic BCs along the x and y directions. Each lattice 
site [x, y, z) with I < x < L^, 1 < y < Ly, 1 < z < 
and lattice constant 1 is occupied by a spin Sx,y,z = il- 
The Hamiltonian of the Ising model with surface fields is 

-J = — Sx^y^zSx',y',z' + Hi Sx,y,l + H^ Sx,y,L^, 

(nn) x,y x,y 

_ (2) 

where J > is the spin-spin interaction constant, H^ = 
H^ J and H^ = H^ J arc the values of the surface 
boundary fields acting on the spins in the bottom and in 
the top layer, respectively. The sum (nn) is taken over all 
nearest-neighbor pairs of sites on the lattice and the sum 
X, y corresponding to the boundary fields is taken over 
the top and the bottom layer. Here we do not consider 
a bulk field. In the following temperatures, the surface 
fields, and energies are measured in units o f J; the inverse 
critical temperature is /3c = 0.2216544(3) [37| . 

For a fixed width Lz and a fixed aspect ratio p = 
Lz/Lx = Lz/Ly of the slab, the thermodynamic state of 
the system is characterized by three parameters: t,H^, 
and H^ . Based on finite-size scaling arguments, for the 
present system Fisher and Nakanishi [38j proposed the 
following convenient scaling variables associated with the 
surface fields: 

hf:^HtLf"^'/^ = «f [L./(^f/^o+)]^°^'/''(3) 
(c^^r^ [Lz/ief/ii^f; 

£i and if correspond to bottom and the top surface, 
respectively. 

Here we study the following three trajectories (see 
Fig.IH): 

(I) hf — oo, an infinitely strong top surface field. 
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(I) hi 

(II) ht 

(III) hi 
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FIG. 1. The parameter space spanned by the scahng vari- 
ables and corresponding to the top and bottom surface 
fields and , respectively (see Eq. ((3])). We investi- 
gate the following paths: (I) hf = 00 (red line) correspond- 
ing to an infinitely strong surface field H^' , (II) /i^ = \hi\ 
(green line), (III) = (blue line). Dashed lines of corre- 
sponding colors denote trajectories, which are equivalent due 
to the exchange symmetry -f^ h^^ . Since Eq. ([2]) does 
not contain a bulk field there is in addition the symmetry 
{h+,h-) o {-h+,~h-). 



(II) hi = \hi I, finite symmetric and antisymmetric sur- 
face fields. 

(Ill) = 0, free boundary conditions at the top 
surface. 

In the simulations, case (I) is realized by fixing all spins 
in tlie top layer 2; = at the value -t-1. For fi- 
nite surface fields it is convenient to replace the sur- 
face field applied at the top (bottom) surface of the 
slab by having this surface layer being linked via mod- 
ified bonds to spins located in an extra layer z = ( 
z = Lz + 1) with the interaction — i/f J2x,y Sx,y,oSx,y,i 
i~H^ Y.x,y Sx,y,L,Sx,y,L,+i); the spius in the extra layer 
z = {z ^ L, + 1) are fixed at the same value -1-1 for 
all x,y. In practice, a surface field, which is finite but 
strong enough to lead to a saturation of the data, can be 
used to mimic the action of an infinite surface field. For 
instance, for > 100 we do not observe any variation 
of our data as function of hf . 

By construction, for all three cases there is only one 
scaling variable associated with the two surface fields. In 
the following we use the notation Hi = Hi and hi = hi . 
This fixes the top surface field H^ in accordance with 

(I) - (III). The plane of parameters {hi = hi, hi) is 
shown in Fig. [TJ We note that the cases (I) and (II) with 
symmetric fields coincide at the point (00,00), the cases 

(II) with symmetric fields and (III) coincide for (0,0), 
and finally for the cases (I) and (III) the point (0, 00) 
coincides with (00, 0). 



We have computed the critical Casimir forces for a 
selection of parameters from sets corresponding to the 
cases (I) , (II) , and (III) which in Fig. [1] are denoted by 
solid lines. Points in the plane (ft,^,/i^) corresponding 
to the cases (I), (II), and (III) which are equivalent due 
to the exchange symmetry hi -H- h'^ arc indicated by 
dashed lines. Due to the absence of a bulk field there is 
also the symmetry (ft,^,/ij~) -O- {—hi,—hi). 

For large areas A, the total free energy 
F{l3,H+,Hi ,Lz,A) of the film of thickness 
can be written as 



F{(3,H+,Hi,Lz,A) 
A 



LJ{l3,H+,Hi,Lz) (4) 



= Lzr^\i3)+r'r{f3,H+,Hi,Lz), 

where /'^""*(/3) is the bulk free energy density at a given 
temperature. The excess free energy f^^ per area con- 
tains two Lz-independent surface contributions in addi- 
tion to the finite-size contribution Hi , Hi , Lz) — 



f^^i/S, Hi,Hi , 00) which vanishes for Lz 



The Lz 



dependence of the latter gives rise to the critical Casimir 
force /c per unit area A and in units of k^T = j3~^: 

fcW,H+,Hi,Lz) = -dr{l3,H+,Hi,Lz)/dLz, (5) 

with the bottom surface field Hi = Hi and the upper 
surface field H^ = {00, |i/]"|,0}, in accordance with (I), 
(II), and (III), respectively. 

For a lattice (lattice quantities are denoted by a "hat" 
"), the derivative in Eq. ([5]) is replaced by a finite differ- 
ence and /c(/3,i) is given by 



fcil3,Hi,L,A) = - 



pAFi/3,Hi,L,A) 



A 



pr'\P). (6) 



with the free energy difference AF{(3, Hi, L, A) — 
F{l3,Hi,L + i,A) - F{l3,Hi,L- i,A). In these three 
expressions the thickness L = Lz — ^ is half-integer, so 
that the rhs is expressed via the free energy difference for 
slabs of integer thicknesses Lz = L+^ and Lz — l= L—^. 
Later on we shall denote by Lz the thickness of the sys- 
tem for which we perform the computations and by the 
half-integer quantity i = — i the variable the critical 
Casimir force depends on. 

From the general theory of finite-size scaling (ssl . Isoj 
and based on renormalization-group analyses [40| we ex- 
pect that in the scaling limit the Casimir force takes the 
universal scaling form 

/c(/3, H+,Hi,L) = L-^i) ((L/6)'/''sign(t), hi,hi 

(7) 

where the scaling function ■d{T = {L/S,Q)^^'^t, hf , hi ) de- 
pends on the spatial dimension D and on the bound- 
ary conditions on the top and bottom surfaces. Here 
^f, = ^±|i|-'' is the bulk correlation length which controls 
the spatial exponential decay of the two-point correlation 
function; are nonunivcrsal amplitudes above (-I-) and 
below (— ) the bulk critical temperature Tc- In the whole 
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A = 




FIG. 2. Bond arrangement for the computation of the free 
energy difference in Eq. (|10p between systems of thickness 
Lz (a) and — 1 (b) (see the main text). The crossover 
Hamiltonian "Her (c) belongs to a system which interpolates 
between those described by the Hamiltonian Ho (a) for the 
system of thickness Lz (for A = 0) and by the Hamiltonian 
Hi (b) for the system of thickness Lz ^ 1 plus a 2D layer of 
area A (for A = 1). 



range of temperatures, we plot the scaling functions us- 
ing the value = 0.501(2) [s^l which is the amplitude of 
the second moment correlation length £^2'^d ; for the Ising 
model 6/$2nd =i 1 for /3 < /3c H. 

Below we shall use the following notations: 
^(^){T,hi) = ?9(r>+ =^ oo,h^), ?9(")(r,/ii) = i?(t,/i+ = 
\h^\,h^), and 'd'-'''\T,hi) = i?(t,/i+ = 0,h^). 



III. NUMERICAL METHOD 

We compute the free energy difference AF{/3, Hi,L, A) 
by using the so-called coupling parameter approach (see, 
e.g., Refs. [Ill and [HI). This is a viable alternative to 
the method used in Ref. [i^, in which a suitable lattice 
stress tensor has been introduced in such a way that its 
ensemble average renders AF. So far, this latter method 
can be implemented only for periodic BC. 

The coupling parameter approach is used in order to 
compute the difference Fi — Fq between free energies 
Fi = — In cxp(— /JHi), i = 0, 1, of models char- 
acterized by two different energies as given by Hamilto- 
nian Hq and "Hi. Such a calculation is successful if the 
configuration space C (i.e., the whole set of spins) is the 
same for both models. In order to implement this ap- 
proach, one introduces an interpolating system with the 
crossover Hamiltonian 



•Hcr(A) - (1 - X)Ho + XHi 



(8) 



As a function of the coupling parameter A £ [0,1], 
'Hcr(A) interpolates between "Hq and Jii as A increases 
from to 1. Accordingly the free energy -Fcr(A) = 
— -g In^^ exp(— /3Hcr(A)) of the crossover system inter- 
polates between Fq and Fi. The sum is taken over all 
spin configurations C of the model, which are the same 
for Fq, Fi, and Fcr- The difference Fi — Fq can triv- 



ially be expressed as Fi — Fq = F^j.{\)dX where F^ 
is the derivative of Fci.(A) with respect to the coupling 
parameter: 



dJ^cr(A) Ec(^i-^o)e 



dX 



Ec e- 



-/3W„(A) 



(AH)cr(A), (9) 



which takes the form of the canonical ensemble aver- 
age (. . .)cr(A) of the energy difference AH = Hi — Hq 
with respect to the crossover Hamiltonian Ha for a given 
value of the coupling parameter A. The energy difference 
(AH)cr(A) can be computed efficiently via MC simula- 
tions of the lattice model characterized by the Hamilto- 
nian Hci- Finally, the difference of free energies is ex- 
pressed as an integral over the mean energy difference 
(see, e.g., Ref. 



Fi-Fq= f (AH)cr(A)dA. 

Jo 



(10) 



According to Eq. ([5]) we are interested in the 
difference AF{f3, Hi, L, A) between the free energies 
F{(3,Hi,L,,A) and F{l3,Hi,L^ - 1,A) (we recall that 
L = Lz — ^ )■ In order to apply the method described 
above for the computation of AF{/3, Hi, L, A) (which 
renders fc (see Eq. (H]))) one identifies the model, the 
Hamiltonian Hq, and the associated configuration space 
C with the corresponding quantities of the model we 
are interested in on the lattice A x Lz (see Fig. EJa)) 
so that Fq{(3,Hi,Lz,A) = F{l3, Hi, L^, A). The final 
system Hi is identified with the slab of area A and 
thickness L^ — 1 plus a two-dimensional layer of size A: 
Fi{p,Hi,Lz,A) = F{I3,Hi,Lz-1,A)+F2dW,A) (see 
Fig. Elb)). Here F2d(/3, A) is the free energy of the iso- 
lated 2D layer of area A. One has to include this 2D layer 
into the consideration in order to maintain the same num- 
ber of spins in the configuration space C for the initial, in- 
termediate, and final models. This layer can be extracted 
from the initial model at any position zq = 1, 2, . . . , 
along the z-direction. It decouples from the rest of the 
lattice upon passing from A = to A = 1, i.e., from 
Fig. m (a) to (b) via (c). The corresponding crossover 
Hamiltonian "Her (A) (but not the result of the integra- 
tion in Ea. ([TU)) ) does depend on the position zq from 
where the 2D layer is extracted. In our simulations we 
use Zq = Lz/2 for even values of L^ and zq = {Lz — l)/2 
for odd values of L^- The explicit expression for the en- 
ergy difference Hi — Hq is 



AH — y ^ {Sx,y,zo-lSx,y,zo + l 

~ Sx,y,zo-lSx,y,zo ~ Sx,y,zoSx,y,zo + l) i 



(11) 



where the three indices {x, y, z) identify a lattice site, the 
sum is taken over all lateral lattice site positions in the 
xy plane, and with a coupling strength J — \ (indicated 
by solid bonds in Figs. [5] (a) and (b); J is absorbed into 
13). The crossover Hamiltonian Her (A) = Hq + XAH 
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is characterized by the couphng constants depicted in 
Fig. IHc). The free energy difference AF (see Eqs. © 
and (jlOp) can be expressed as 



{AH)cr{X)dX + F2DW,A) (12) 



AF{p,Hi,L,A) 



where the integral is taken for fixed values of f3 and Hi. 
Note that although AH is independent of Hi, the de- 
pendence of AF on Hi enters via the statistical weight 
- cxp(-;3-Hcr). The free energy F2u(/3,A) of the 2D 
layer can be computed from the analytical expressions 
given in Ref. j43j . 

Once AF{(3, Hi , L, A) has been computed, one has still 
to subtract /'^""^(/3) from it (see Eq. ^) in order to 
obtain the Casimir force for a slab of assigned thickness 
L = Lz — 1/2. We determine the bulk free energy by using 
the temperature integration method [l^, [13, Ejl applied 
to a cubical system of size Lcubc with periodic boundary 
conditions. For such a system the free energy per site (in 
units of ksT) can be written as 



/?/(/?, icubc) = -ln(2) + 



"^cubc 



(£;(/?', icubc))d/3', 



(13) 

where (£'(/?', Lcube)) is the averaged internal energy of 
the system at the inverse temperature and for the 
size Lcube; — lii(2) is the free energy in units of ksT 
and per site at /3 = 0. For a cube, in the limit 
Lcube — ^ OO the finite-size dependence of the free en- 
ergy density for a cube is predicted [sl] to scale with 
Lcube as /3/(/3, Leube) " /?/^"''^(/3) (X L-^^^^. Therefore the 
bulk free energy per spin follows as the limit /3/''""^(/3) = 



lim 



Pf{l3, Lcubc) • At the critical point one has 



/3c/(/?c,icubc)=^/?c/'^"'^/3c 



UnL 



-3 



O^cubC 



(14) 



with the universal finite-size scaling amplitude Uq — 
-0.657(3) (see Ref. 

In order to determine the universal scaling function 
of the critical Casimir force we perform the following 
steps (details are given below). For each temperature 
we compute the averaged internal energy {E{f3, Lcuhc)) 
for a cube with periodic boundary conditions by using 
a histogram reweighting MC method. Then we carry 
out a numerical integration in order to obtain an esti- 
mate for the bulk free energy /3/'^"''^(/3) in accordance 
with Eqs. ([T3| and (fT4)) . For the slab geometry A x L^, 
at the inverse temperature /3, and for a fixed boundary 
field Hi, we compute the ensemble averages (A'H)cr(A) 
via MC simulations for N\ = 21 different values of 
k = 0, . . . ,N\-i. Based on these N\ val- 



ues we carry out the numerical integration in Eq. (|12p 
and use an analytical expression for F2d{(3,A) in order 
to obtain AF{f3, Hi, L, A). Combining the results for the 
bulk free energy density /3/'^""*^(/3) and for the free energy 



difference and by using Eq. ([6]) we obtain a numerical es- 
timate for the critical Casimir force fc{j3. Hi, L, A). In 
order to obtain the corresponding scaling function ■& we 
perform computations for various values of L, A, the in- 
verse temperature /?, and boundary fields Hi. The scal- 
ing function 'd in Eq. ([7]) is retrieved from the numerical 
data for fc by taking into account finite-size corrections 
as described in the following section. 

For determining the bulk free energy density the his- 
togram reweighting method has been used as follows 
psL l46j . The computation of the energy distribution 
P{E,/3i) has been performed for a choice of 256 points 
Pi g [0, 0.3] for a cubic system of size Lcubo = 128. For 
the numerical simulation we have employed the hybrid 
MC method, which is a suitable mixture of Wolff and 
Metropolis algorithms [i^. For thermalization 4 x 10'^ 
hybrid MC steps have been used. The averaging has been 
performed over 10^ hybrid MC steps which have been 
split into 10 series for the evaluation of statistical errors. 
Therefore, for every value Pi actually ten histograms 
(each consisting of 10^ MC steps) have been computed. 
According to the histogram reweighting method one can 
obtain an estimate for (E) at an inverse temperature 
based on the histogram P{E,j3i) for the inverse temper- 
ature ft Hi,!!!: 



j:EP{E,(3,)e~^(f''-^') 

_E 

j:F{E,P,)e~EiP'-0,) ■ 

E 



(15) 



For every /?' e [/3i,/3i+i] we define the interpolated inter- 
nal energy 



ft 



ft 



Pi+l — Pi 

(16) 

We have checked that for the same inverse temperature 
/?' the difference between the estimates {E)p^{[i') and 
(E) p.^-^{P'), which use histograms for two neighboring 
points ft and ft+i, is substantially less than the statis- 
tical inaccuracy of our simulation data. The statistical 
inaccuracy has been determined canonically over 10 se- 
ries of histograms. In the next step, in accordance with 
Eq. we obtain the free energy /3/(/3, Lcube) by inte- 
grating numerically the interpolated internal energy. For 
the intergration we employ the trapezoidal rule with a 
large (> 10^) number of points, so that the inaccuracy 
of the numerical integration is less then 10~^. We es- 
timate that at the bulk critical temperature Pc the sta- 
tistical error APcfiPc, Lcube) for the free energy deter- 
mined from 10 series is about 4 x 10^*". In the follow- 
ing, we neglect the finite-size correction of the bulk free 
energy and take /3/*'""^(/3) :^ /3/(/3, Lcube = 128) (com- 
pare Eq. ^^). This is justified, because for the maximal 
value L = 19.5 used in our simulation the finite-size cor- 
rection to 1? due to the finite system size Lcube = 128, 
i.e., 0. 657 (L/Lcube)^ — 0.0023, is of the same order 
as the statistical error stemming from the contributon 
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L^ApjiP,, L,ube), i.e., (19.5)3 x 4 x lO'^ ~ 0.00297 (see 
Eqs. (El) and 0). 
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FIG. 3. MC data for case (I): (a) Casimir force fc as a 
function of the inverse temperature j3; /3c = 0.2216544(3) (b) 
rescaled Casimir force fc as a function of the scahng vari- 
able r = [L/^^y^'^iT - n)/Tc = {L/i, + Y'"t. The data 
correspond to hi = -100,0,100 and L = 9.5,14.5,19.5. In 
(a) and (b) the data for (/ii, L) = (100, 9.5) and (100, 14.5) as 
well as (0, 14.5) and (0, 19.5) can be barely distinguished. 

For the computation of the free energy 
Af'(/3,i?i,i, A) in Eq. (O we use slabs of thick- 
nesses Lz = 10, 15, 20, so that L = 9.5, 14.5, 19.5, with 
an aspect ratio equal to 6: = Ly = A = 36i^. 
In order to compute the average (A'H)cr(A) we again 
use the hybrid MC method with a mixture of Wolff and 
Metropolis algorithms. Each hybrid MC step consists of 
a flip of a Wolff cluster according to the Wolff algorithm, 
followed first by 3A attempts to flip an arbitrary spin 
and then by 3A attempts to flip a spin s^.y.z with 
z S {zq — I, zq, zq + 1}. These attempts are accepted 
according to the Metropolis rate H^. We use 2.5 x 10^ 
MC steps for thermalization. For the computation of the 
thermal average we use 5 x 10^ MC steps split into 10 
series. For each series, using Simpson's rule we perform 
a numerical integration over N\ = 21 points for fixed 
values of the inverse temperature /3, the surface field 



Hi, and the width L of the slab. Having computed 
the free energy difference AF{(3, Hi, L, A), for each 
series we finally combine the results for the bulk free 
energy /3/'^^'^(/3) with the corresponding ones for the 
free energy difference AF{(3, Hi, L, A) and determine 
the numerical inaccuracy. 

In Fig. EJa) we plot the Casimir force fc as a function 
of /3 for the three values hi = —100, 0, 100 of the bottom 
surface scaling field. In Fig.[3][b) we plot the rescaled val- 
ues of the Casimir force L^fc as a function of the tem- 
perature scaling variable r = {L/^+^/^iT - Tc)/Tc = 
{L/^^^Y/'^t for case (I). The visible absence of the ex- 
pected data collapse is due to finite size corrections to 
scaling, which will be discussed in the following section. 



IV. FINITE SIZE CORRECTIONS TO SCALING 

Finite-size scaling is known to be valid asymptotically 
for finite but large lattices and small values of t, i.e., for 
a large bulk correlation length here large means rela- 
tive to the lattice constant [s^ . Outside the asymptotic 
regime corrections to the leading (universal) scaling be- 
havior become relevant. These non-universal corrections 
affect both the scaling variables and the scaling func- 
tions and depend on the details of the model as well as 
on the geometry and the boundary conditions [13, 13 ■ 
Renormalization-group analyses reveal that there is a 
whole variety of sources for corrections to scaling which 
arise from bulk, surface, and finite-size effects [39| . 

For the finite and rather limited sizes of the lattices 
which wc investigate in our MC simulations, it is neces- 
sary to take corrections to scaling into account in order 
to obtain data collapse and thus allowing us to infer the 
leading universal scaling function pi]| - [l2| . 

In the present study, the following quantities are ex- 
pected to acquire corrections to scaling: 

• the amplitude of the scaling function d ~ fc 

• the surface field scaling variable hi 



• the temperature scaling variable t ~ t (i/Co^) 



1/^ 



In our previous MC simulations aimed at obtaining criti- 
cal Casimir forces for Ising films with a variety of univer- 
sal boundary conditions, such as (-I-, -I-), (-I-, — ), or (O, O) 
BC [13, [ill, corrections to scaling were taken into ac- 
count by using various ansatze. The choice for a particu- 
lar form of corrections to scaling was guided by achieving 
the best data collapse or the best fits used in our com- 
putations. For example, for the amplitude of the scaling 

function we adopted the expression fc = ^~^ ^(iX%^--^) ^- 
Various variants for this form of corrections to scaling 
were considered; we used {gi ^ 0, 52), (ffi, 52 7^ 0), or 
[gi 7^ 0, (72 7^ 0). They all lead to a satisfactory data col- 
lapse, but the inferred amplitude of the scaling function 
of the critical Casimir force depends sensitively on the 
particular ansatz. 
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FIG. 4. Results at the critical point P — Pc for case (I) (see 
Fig. [T]) as a function of boundary fields: (a) Casimir force 
multiplied by L'^ as a function of the scaling variable hi = 
HiL^^ without finite size corrections; (b) Casimir force 
scaling function i5'^'(f = 0, h^) with taking corrections L + 5 
into account as a function of the corrected scaling variable 
hi — Hi{L + 5)"^! in units of the lattice spacing 5 = 0.65 
at the bulk critical point /3c. For hi = one finds the fixed- 
point (O, -f ). 



In Rcfs. [3J, |4J] still another type of finite-size eor- 
rection is employed. It amounts to introdueing an ef- 
fective width L + S oi the slab so that, e.g., the ampli- 
tude of the scaling function of the critical Casimir force 
scales as fc = {L + S)~^i). Using this type of finite- 
size correction may be justified as follows. As mentioned 
in Sec. |T] surfaces subjected to the action of a surface 
field asymptotically belong to the surface universality 
class of the normal transition which corresponds to (+) 
or (— ) fixed-point boundary conditions in the sense of 
renormalization-group theory [l^, [2l|. (The (-I-) and (— ) 
boundary conditions are realized as the limits of the scal- 
ing field /i^'' — >■ -l-oo and — oo, respectively.) For such 
boundary conditions, on the coarsed-grained scale the or- 
der parameter varies as \(j){z — 0)| oc z~^l^ for small nor- 
mal distances from the surface (but still large on molecu- 
lar scales) 0, . Within a certain range of small z val- 



ues such a divergent behavior is expected to hold also for 
a finite but sufficiently strong surface field. In Ising lattice 
models with boundary conditions corresponding to (4-) 
or (— ) fixed-point BCs, the order parameter does not di- 
verge at the surface but saturates there at the value -1-1 or 
— 1. Changing the width of the slab from LtoL-Vb with 
a nonuniversal length i5 = z^ -f z^ such that the order 

parameter profile behaves as |0(z — J- 0)| '--^ (z -I- z^cl)~^^^ 
[20l . [sot upon approaching the wall z, turns out to be 
an effective means to take into account corrections to 
the leading critical behavior [sot: Zex plays the role of 
an extrapolation length [^O, l2ll|. Similarly, the effects of 
a physical wall with a finite surface field (which implies 
6^ 7^ ( sec Eq. H]))) on the order parameter are equiva- 
lent to those of a fictitious wall with strong surface fields 
(which means if' = 0) displaced by a distance — Zox 
from the physical wall. One can determine the length b 
by analyzing the spatial variation of the order parame- 
ter profile, as it was done for the Blume-Capcl model in 
Ref. [l^l- Here, we assume that the equivalence described 
above carries over to critical Casimir forces such that we 
can determine the effective width L -f (5 of the slab by de- 
manding the best data collapse. We apply this method 
also in the crossover regime, i.e., for sufficiently weak sur- 
face fields for which upon approaching the critical point 
one effectively observes a crossover to the boundary con- 
dition corresponding to the ordinary transition (O) fixed 
point. As discussed in the introduction, the order pa- 
rameter profiles in a film with weak surface fields deviate 
strongly from the fixed-point universal behavior. Accord- 
ingly, we expect that within this range of surface fields 
the aforementioned type of correction does not satisfac- 
torily capture the actually corrections to scaling. 

In Fig. mja) we plot the rescaled critical Casimir force 
^^]c for case (1). It is evaluated at the critical point 
/3c and presented as a function of h\ without finite-size 
corrections taken into account. Apparently the data for 
the rescaled force do not coincide for various values of 
L = 9.5, 14.5, 19.5, 24.5. In order to obtain the expected 
data collapse we apply the following finite-size corrections 
(here and in the following we denote scaling variables 
with finite size corrections by a tilde: f , /ii): 



with 



f ^r[(L + <5)/^+]'/'' [1 + 5„(L + ^)- 
and (see Eq. ^) 

^1 =i7i(L + <5)^""/^ 



(17) 



(18) 



(19) 



where cj = 0.84(4) is the leading bulk correction-to- 
scaling exponent [2J]; the length b and the coefficient 
remain to be determined. 

The value of the length i5 is obtained from the data for 
the critical Casimir force at the critical point (these data 
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are presented in Fig. lUJa)). By using Eqs. ([TT)) and ^T9\\ 
and by implementing the fitting procedure within the in- 
terval hi e [—15, 15] with S being the only fit parameter, 
we obtain the value 6 = 0.65 which minimizes deviations 
between data for different values of L. Including error 
bars we find S ~ 0.65(2) for various intervals of hi or 
S = 0.60(5) for different sets of L. The final result for 
the scaling function '!9(0, /ii) with corrections to scaling 
corresponding to S = 0.65 is shown in Fig. IH^b). For 
large absolute values of hi we reproduce the data from 
Refs. [l3i[ll| for the critical point with (— , +) and (+, +) 
BCs. For small values of \hi\ we observe the crossover 
between these two regimes. 



various -dk: -dexpected = (1/^) Z^fcLi ^^fe- Finally, for ev- 
ery Lfe and for a given value of 5 we compute the sum of 
squares {5) of the deviation of the aforementioned scal- 
ing functions from the expected scaling function. Finally, 
we determine as the value of 5 the one which minimizes 
. That value provides the best data collapse of the 
data for different L. 

Applying the same procedure for case (II) and case 
(III) we obtain the values 5 = 0.6 and 5 = 1.4, respec- 
tively. However, in order to be consistent (as mentioned 
earlier, for some values of the surface field different cases 
coincide) we use the common value 5 = 0.65 for all cases. 
In Figs. [5] and ini we present our results without (a) and 
with (b) finite-size corrections for case (II) and case (III), 
respectively. 
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FIG. 5. Same as Fig. g] for the case (II) (see Fig. [T]). For 

fti = one finds the fixed point (O, O). 

The procedure which we used in order to obtain the 
best fit for the value of the length 5 is described in detail 
in the appendix of Ref. [llj. One of the difficulties in 
finding the optimal data collapse is that the fitting func- 
tion itself, i.e., the scaling function of the critical Casimir 
force, is not known. For the initial guess for the value of 
the length 5 we infer the scaling function from the corre- 
sponding data for /c, one function 'd^ for each value 
(fc = 1, • • • , A^) used (see Eqs. (d?]), (HH]), and dH])). Wc 
define an expected scaling function as the average of the 
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FIG. 6. Same as Fig. g] for the case (III) (see Fig. [T]). For 
/ii = one finds the fixed point (O, O). 



Knowing the finite-size corrections of the surface field 
scaling variable we can carry out numerical simulations 
for various values oi hi] for each value oi hi we can ex- 
tract information about the coefficient g,^ using the same 
procedure as for the determination of the length 5. 
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FIG. 7. The scaling function of the critical Casimir 
force for case (I), i.e., — oo as a function of the tem- 
perature scaling variable f (see Eq. (|18p ) for various bot- 
tom boundary fields corresponding to certain values of the 
surface field scaling variable hi = h~[ (see Eq. (|19|l '): (a) 
large amplitudes of the surface field (from top to bottom): 
hi = -100,-8,-4,-2,0,2,6,100; (b) small amplitudes of 
the surface field for which a crossover from repulsive to at- 
tractive forces as function of f is observed (from top to bot- 
tom): h\ = —1, 0, 0.5, 1, 1.5, 2. For each color corresponds 
to L = 9.5, ■ to L = 14.5, and A to L = 19.5. 



V. RESULTS 



FIG. 8. Same as Fig. [7] for case (II) (see Fig. [l]). In (b) 
for small amplitudes of the surface field [from top to bottom: 
hi = h^ — —2, —1, 0, 0.5, 1, 1.5] there occurs a crossover from 
attractive to repulsive forces upon increasing f. 



TABLE I. Values of hi with the corresponding values of g^i as 
obtained from the fitting procedure for case (I), i.e., hf = oo. 



hi 


-100 


-8 


-4 


-2 


-1 







-0.56(2) 


-0.05(2) 


-1.19(2) 


1.04(2) 


1.60(4) 


2.3(2) 


hi 


0.5 


1 


1.5 


2 


6 


100 




0.68(12) 


0.72(15) 


-0.05(5) 


-0.145(3) 


-0.47(2) 


-0.94(2) 



Here we present the critical Casimir force scaling func- 
tion determined for various values of hi as a function of 
f. The set of values used for is given in Tables HI HIl 
and mil For each of the three values L = 9.5, 14.5, 19.5 
of the slab thickness we infer the values of the surface 
fields Hi which correspond to the pair {hi,L) according 
to Eq. ([TO)) with S = 0.65. Next, for each pair {hi,L) the 
critical Casimir force /c(/?, Hi,L) has been computed for 
various inverse temperatures /?. Finally, for each value of 
hi we apply the fitting procedure described above in or- 
der to determine by using Eq. ([T5|) . Our results for 
guj are given in Tables HI HIl and lllll for the cases (I), (II), 
and (III), respectively. 



TABLE II. Same as TableUfor case (II), i.e., h+ = \hi 



hi 


-100 


-8 


-4 


-2 


-1 







-0.58(2) 


-0.14(2) 


-0.25(2) 


0.18(2) 


0.10(3) 


0.70(3) 


hi 


0.5 


1 


1.5 


2 


6 


100 


goj 


2.05(10) 


0.51(10) 


0.97(13) 


1.30(5) 


-0.06(2) 


-0.95(2) 



TABLE III. Same as TableUfor case (III), i.e., h+ = 0. 



hi 


100 


6 


2 


1 


0.5 







2.8(1) 


1.6(1) 


-0.12(5) 


1.51(5) 


1.87(4) 


1.30(3) 
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FIG. 9. Same as Fig. [7] for case (III) (/ii = h~ see Fig. [TJ. 



For a selection of values of the surface field scafing vari- 
able hi in Fig. ^a) we present our results for the critical 
Casimir force scaling function i?'-^^(f , /ii) corresponding 
to case (I). We find that de facto for hi = ±100 the 
scaling limit of infinitely strong surface fields has been 
reached and the scaling function of the critical Casimir 
force corresponds to the fixed-point BCs (-|-,±). In 
Fig. [TJb) we present our results for small values of \hi\, 
hi G [—1,2], for which we observe a transition from a re- 
pulsive to an attractive force upon increasing the temper- 
ature scaling field f. In these instances in which a change 
in sign is observed in Fig. ^h) , the length scales associ- 
ated with the surface fields are £f = (see Eq. ([T]) for 
= oo) on the top and Kr/Co"! = 0.44L^ (see Eqs. ^ 
and © for hi = 1 and c = 0.5/^^) and lif/^^l ~ 1.13L^ 
(see Eqs. ^ and ^ for hi = 0.5 and c = 0.5/^^) on the 
bottom. (In D ~ 3, c ~ 1/a ii the coupling constant in 
the surface row is unchanged relative to the one in the 
bulk 0,[2l|; with 5+ ~ 0.501 this implies ~ 0.501.) 

In Fig. Ufa) we show data for the critical Casimir force 
scaling function (f, ft-i) corresponding to case (II). 
For hi = 100 and hi = —100 we recover the scal- 
ing limit corresponding to (++) and (H — ) fixed-point 
BCs, respectively. We note that the change in sign of 
the critical Casimir force upon varying the temperature 
occurs only for opposing surface fields. As before this 
change of sign is observed for weak surface fields for which 
\ei/^+\^L,,i.e.JoT hi = -2,-1. 

Finally, in Fig. [S] the data for case (III) are presented. 
This case contains in particular two fixed-point BCs: for 
the value hi = 100 we observe a universal behavior of the 
scaling function of the critical Casimir force correspond- 
ing to (O, fixed-point BCs, whereas for ft-i = we find 
the (O, O) fixed-point universal behavior of i?(f , hi). As 
in the other cases, the crossover from attraction to re- 
pulsion can be achieved by increasing the temperature 
scaling variable f, provided that the surface fields are 
sufficiently weak so that l^'^/Co'l ~ L^, i.e., for hi = 2 
and 1. 



VI. SUMMARY AND CONCLUSIONS 

For variable surface fields we have determined via MC 
simulations the universal scaling functions •& of critical 
Casimir forces for 3-D Ising slabs describing the crossover 
from the ordinary to the normal surface universality class 
(Figs. [71 [HI and [HI). This amounts to investigate the scal- 
ing functions i?(r, ft,^, /ij~) (see Eq. ^) for finite values 
of the surface fields. We have computed the lattice scal- 
ing functions i?(f, ft,i) along three different paths in the 
parameter space {h'^,h^) (see Fig. [1]): ??(^'(f, /ii) cor- 



responding to /ij" 



oo, 



(f, /ii) corresponding to 



h^ = \hj^ I, and i5'^^^l^(f , /ii) corresponding to h'^ = 0. 
Due to the fact that on the lattice the derivative in Eq. ([5]) 
is replaced by a finite difference, the scaling function 
as function of the corrected scaling variables {f,hi) esti- 
mates the leading behavior of d as function of (r, h'^ ,h^); 
alternative definitions of the lattice derivative give rise to 
distinct corrections for both the scaling function and the 
scaling variables. We have focused on cases in which 
upon variation of the temperature a crossover from at- 
tractive to repulsive critical Casimir force is observed. 
Such a behavior is particularly interesting in view of po- 
tential application, e.g., for colloidal suspensions. We 
have found that a change of sign of the critical Casimir 
force as a result of a minute change in temperature oc- 
curs only in systems with strongly asymmetrical surfaces, 
i.e., in cases in which the two surface fields differ signif- 
icantly in magnitude. For this phenomenon to occur at 
least one of the surface fields has to be weak enough such 
that the length scale £i associated with the surface field 
Hi (Eq. (|T])) is comparable with the width L of the slab 
(see Figs. [71[b) and M corresponding to fixed h'^ = oo 
and h^ = 0, respectively, and a variable second surface 
scaling field hi). We note that for such large values of 
£i the order parameter profiles near a single wall differ 
significantly from the ones corresponding to strong sur- 
face fields which belong to the surface universality class 
of the normal transition. If both surface fields are weak 
and have the same magnitude they must have opposite 
signs in order to produce a change of sign of the critical 
Casimir force (see Fig. IM^b)). The change from attrac- 
tion to repulsion (i.e., a zero of "d) can occur either below 
the bulk critical temperature, as for the cases in which 
one of the surfaces is subjected to the (O) fixed-point 
BC or for weak opposing surface fields (see Figs. [HI and 
[SI respectively), or above Tc, as for the (-I-) fixed-point 
BC (see Fig. [Tl^b)). In all cases the change of sign takes 
place rather close to the critical point. 

Corrections to scaling have had to be taken into ac- 
count in order to obtain data collapse which allowed us 
to infer the universal scaling functions (see Fig. [3]). The 
introduction of an effective width L + S oi the slab turned 
out to be a very useful way of implementing corrections to 
scaling, provided the surface fields are not too weak. The 
value of the length S has been obtained from the data for 
the critical Casimir force at the critical point (sec Figs.[3| 
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and |6] corresponding to fixed values h'^ = oo and hf = 
for the top surface, respectively, and a variable surface 
scaling field for the bottom surface, and Fig. [5] corre- 
sponding to the surface fields for the two surfaces to be 
of the same magnitude; compare Fig. [1]). 

The present results close an important gap in the 
knowledge of the Casimir scaling function for the 3D 
Ising universality class. The theoretical results for vari- 
able surface fields have been available in D = 2 (from ex- 
act calculations in Ising strips [s^l ) and in D = 4 (from a 
field-theoretic approach [s^] ) ■ The MC simulation results 
in Ref. [s^l have been obtained for the 3D Blume-Capel 
model which is an extension of the Ising model studied 
here. They provide critical Casimir forces as function of 
Pc — P for certain values of the surface fields, which in the 
parameter space shown in Fig. [T] correspond to path (/) 
with > 0. Because the choice of the surface fields for 
the presented data is different from ours, we cannot make 



a direct quantitative comparison (except for the case of 
{h^ = oo,h^ = 0) for which the data agree). However, 
there is a qualitative agreement with our findings; for cer- 
tain choices of the surface fields the critical Casimir force 
changes sign as function of temperature. This agreement 
provides further evidence for the universal character of 
critical Casimir forces. 

Our data for the critical Casimir scaling function have 
the crucial advantage over the results in D = 2 and 4 that 
they can be directly compared with possible experimen- 
tal data. Interestingly, in all spatial dimensions studied 
the crossover behavior of the scaling function of the crit- 
ical Casimir force as a function of the temperature scal- 
ing variable is qualitatively the same. The robustness of 
this observation indicates that an experimental observa- 
tion of the change of sign of the critical Casimir force 
with temperature is possible, provided that the chemical 
properties of the confining surfaces are carefuly chosen. 



[1] M. E. Fisher and P. G. de Gennes, C. R. Acad. Sci. Paris 
Ser. B 287, 207 (1978). 

[2] M. Krech, Casimir Effect in Critical Systems (World Sci- 
entific, Singapore, 1994); J. Phys.: Condens. Matter 11, 
R391 (1999). 

[3] J. G. Brankov, D. M. Dantchev, and N. S. Tonchev, The 
Theory of Critical Phenomena in Finite-Size Systems - 
Scaling and Quantum Effects (World Scientific, Singa- 
pore, 2000). 

[4] for a recent review see A. Gambassi, J. Phys.: Conf. Ser. 
161, 012037 (2009). 

[5] U. Nellen, J. Dietrich, L. Helden, S. Chodankar, K. 
Nygard, J. F. van der Veen, and C. Bechinger, Soft Mat- 
ter 7, 5360 (2011). 

[6] R. Garcia and M. H. W. Chan, Phys. Rev. Lett. 83, 1187 
(1999); A. Ganshin, S. Scheidemantel, R. Garcia, and M. 
H. W. Chan, Phys. Rev. Lett. 97, 075301 (2006). 

[7] R. Garcia and M. H. W. Chan, Phys. Rev. Lett. 88, 
086101 (2002). 

[8] M. Fukuto, Y. F. Yano, and P. S. Pershan, Phys. Rev. 

Lett. 94, 135702 (2005). 
[9] S. Rafai, D. Bonn, and J. Meunier, Physica A 386, 31 

(2007). 

[10] O. Vasilyev, A. Gambassi, A. Maciolek, and S. Dietrich, 

EPL 80, 60009 (2007). 
[11] O. Vasilyev, A. Gambassi, A. Maciolek, and S. Dietrich, 

Phys. Rev. E 79, 041142 (2009). 
[12] A. Hucht, Phys. Rev. Lett. 99, 185301 (2007). 
[13] R. Zandi, J. Rudnick, and M. Kardar, Phys. Rev. Lett. 

93, 155302 (2004); R. Zandi, A. Shackell, J. Rudnick, M. 

Kardar, and L. P. Chayes, Phys. Rev. E 76, 030601 (R) 

(2007) . 

[14] A. Maciolek, A. Gambassi, and S. Dietrich, Phys. Rev. 

E 76, 031124 (2007). 
[15] A. Maciolek, and S. Dietrich, Europhys. Lett. 74, 22 

(2006). 

[16] M. Krech, Phys. Rev. E 56, 1642 (1997). 

[17] Z. Borjan and P. J. Upton, Phys. Rev. Lett. 101, 125702 

(2008) . 

[18] C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, and C. 



Bechinger, Nature 451, 172 (2008). 
[19] A. Gambassi, A. Maciolek, C. Hertlein, U. Nellen, L. 
Helden, C. Bechinger, and S. Dietrich, Phys. Rev. E 80, 
061143 (2009). 

[20] K. Binder, in Phase Transitions and Critical Phenom- 
ena, edited by C. Domb and J. L. Lebowitz (Academic, 
London, 1983), Vol. 8, p. 1. 

[21] H. W. Diehl, in Phase Transitions and Critical Phenom- 
ena, edited by C. Domb and J. L. Lebowitz (Academic, 
London, 1986), Vol. 10, p. 76. 

[22] H. W. Diehl, Int. J. Mod. Phys. B 11, 3503 (1997). 

[23] R. Guida and J. Zinn Justin, J. Phys. A: Math. Gen. 31, 
8103 (1998). 

[24] A. Pelissetto and E. Vicari, Phys. Rep. 368, 549 (2002). 
[25] U. Ritschel and P. Czerner, Phys. Rev. Lett. 77, 3645 

(1996); P. Czerner and U. Ritschel, Int. J. Mod. Phys. B 

11, 2075 (1997); P. Czerner and U. Ritschel, Physica A 

237, 240 (1997). 
[26] A. Maciolek, A. Ciach, and J. Stecki, J. Chem. Phys. 

108, 5913 (1998). 
[27] A. Ciach and U. Ritschel, Nucl. Phys. B 489, 653 (1997). 
[28] R. Z. Bariev, Theo. Math. Phys. 77, 1090 (1988). 
[29] A. Maciolek, A. Ciach, and A. Drzewiriski, Phys. Rev. E 

60, 2887 (1999). 
[30] A. Maciolek, R. Evans, and N. B. Wilding, J. Chem. 

Phys. 119, 8663 (2003). 
[31] A. Maciolek, A. Drzewihski, and P. Bryk, J. Chem. Phys. 

120, 1925 (2004). 
[32] T. F. Mohry, A. Maciolek, and S. Dietrich, Phys. Rev. E 

81, 061117 (2010). 
[33] D. B. Abraham and A. Maciolek, Phys. Rev. Lett. 105, 

055701 (2010). 
[34] M. Hasenbusch, Phys. Rev B 83, 134425 (2011). 
[35] F. Parisen Toldin and S. Dietrich, J. Stat. Mech., P11003 

(2010). 

[36] U. Nellen, L. Helden, and C. Bechinger, EPL 88, 26001 
(2009). 

[37] C. Ruge, P. Zhu, and F. Wagner, Physica A 209, 431 
(1994). 

[38] M. E. Fisher and H. Nakanishi, J. Chem. Phys. 75, 5857 



13 



(1981). 

[39] M. N. Barber, in Phase Transitions and Critical Phenom- 
ena, edited by C. Domb and J. L. Lebowitz (Academic, 
New York, 1983), Vol. 8, p. 149; V. Privman, in Finite 
Size Scaling and Numerical Simulation of Statistical Sys- 
tems, edited by V. Privman (World Scientific, Singapore, 
1990), p. 1. 

[40] M. Krech and S. Dietrich, Phys. Rev. Lett. 66, 345 
(1991); Phys. Rev. A 46, 1886 (1992); Phys. Rev. A 46, 
1922 (1992). 

[41] K. K. Mon, Phys. Rev. B 39, 467 (1989); K. K. Mon and 

K. Binder, Phys. Rev. B 42, 675 (1990). 
[42] D. Dantchev and M. Krech, Phys. Rev. E 69, 046119 

(2004). 



[43] B. Kaufman, Phys. Rev. 76, 1232 (1949). 

[44] M. Hasenbusch, Phys. Rev. B 82, 104425 (2010). 

[45] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 

63, 1195 (1989). 
[46] D. P. Landau and K. Binder, A Guide to Monte Carlo 

Simulations in Statistical Physics (Cambridge University 

Press, London, 2005), p. 155. 
[47] V. Privman and M. E. Fisher, J. Phys. A 16, L295 (1983). 
[48] J. M. Luck, Phys. Rev. B 31, 3069 (1985). 
[49] S. Leibler and L. Peliti, J. Phys. C: Solid State Phys. 

14, L403 (1982); E. Brezin and S. Leibler, Phys. Rev. B 

27, 594 (1983). 

[50] M. Smock, H. W. Diehl, and D. P. Landau, Ber. Bunsen- 
ges. Phys. Chem. 86, 486 (1994). 



